function s_mark(N)

  close all
  

  %time_steps = [0:N-1]./(N-1)*1.0;


  function hsp = hsourcep(z)
    if z < 0.001
      hsp = 6+(-8).*z+(-7/3).*z.^2+17.*z.^3+(-33/2).*z.^4+(-6).*z.^5+(421/14).*z.^6+(-181/7).*z.^7+(-94/9).*z.^8+(4672/105).*z.^9+(-2507/70).*z.^10+(-139/9).*z.^11+(16360/273).*z.^12+(-21036/455).*z.^13+(-1879/90).*z.^14+(27723/364).*z.^15+(-207577/3640).*z.^16+(-80/3).*z.^17+(183899/1976).*z.^18+(-2356147/34580).*z.^19+(-2293/70).*z.^20;
    else
      hsp = (1/3).*z.^(-1).*(1+z+z.^2).^(-2).*(3.^(1/2).*pi.*(2+z)+3.*z.*(2+z.*(2+z))+(-2).*3.^(1/2).*(2+z).*atan(3.^(-1/2).*z.^(-1).*(2+z))+3.*(2+z).*log(1+z+z.^2));
    end
  end
  function ksp = ksourcep(z)
    if (z < 0.001)
      ksp = (14/3)+(25/3).*z+(-22/3).*z.^2+(-7/3).*z.^3+(221/21).*z.^4+(-325/42).*z.^5+(-244/63).*z.^6+(7741/630).*z.^7+(-5167/630).*z.^8+(-217/45).*z.^9+(12323/910).*z.^10+(-23473/2730).*z.^11+(-357/65).*z.^12+(39601/2730).*z.^13+(-48787/5460).*z.^14+(-1081/180).*z.^15+(475747/31122).*z.^16+(-5742899/622440).*z.^17+(-76843/11970).*z.^18+(7276627/456456).*z.^19+(-21640043/2282280).*z.^20;
    else
      ksp = (1/3).*z.^(-3).*(1+z+z.^2).^(-1).*(3.*z.*(2+z).*(4+z.*(4+3.*z))+3.^(1/2).*pi.*((-4)+z.*(1+z).*((-6)+z.*(2+z.*(2+z))))+(-2).*3.^(1/2).*((-4)+z.*(1+z).*((-6)+z.*(2+z.*(2+z)))).*atan(3.^(-1/2).*z.^(-1).*(2+z))+3.*((-4)+z.*(1+z).*((-6)+z.*(2+z.*(2+z)))).*log(1+z+z.^2));
    end
  end


  function resMh = Mh(z)
    resMh = zeros(2,2);
    resMh(1,1) = 1;
    resMh(2,1) = 4;
    resMh(2,2) = z;
  end
  function resh = deqh(z, u)
    h   = u(1);
    hz  = u(2);
    resh = zeros(2,1);
    resh(1) = hz;
    resh(2) = hsourcep(z);
  end

  function resMk = Mk(z)
    resMk = zeros(1,1);
    resMk(1,1) = 2;
  end 
  function resk = deqk(z, u)
    k  = u(1);

    resk = zeros(1,1);
    resk(1) = ksourcep(z);
  end


  z = chebpoints(N, 0, 2);

  
  na2t2 = sa_mark(N, 2);
  na2t4 = sa_mark(N, 4);
  na2t5 = sa_mark(N, 5);


  % ******* h k ******

  reltol = 1e-8;

  uin = [0; 6./4];
  opts = odeset('RelTol', reltol, 'AbsTol', 1e-10, 'MStateDependence', 'none', 'Mass', @Mh);
  [tout, uout] = ode15s(@deqh, z, uin, opts);
  nh2s7 = uout(:, 1);
  nh2s7 = 1/2 + z'.*nh2s7; % with z^2 factored out


  uin = [0]; % initial condition: k = 0, controls the 1/z piece
  opts = odeset('RelTol', reltol, 'AbsTol', 1e-10, 'MStateDependence', 'none', 'Mass', @Mk);
  [tout, uout] = ode15s(@deqk, z, uin, opts);
  nk2s7 = uout(:, 1);
  nk2s7 = 2 - 2.*log(z').*z' + z'.*nk2s7; % nk2s7mark with 1/z^2 factored out
  nk2s7(1) = 2; % get rid of the NaN

  save('nfmark.mat', 'nh2s7', 'nk2s7', 'na2t2', 'na2t4', 'na2t5');

  figure
  plot(z, na2t2, z, na2t4, z, na2t5, z, nh2s7, z, nk2s7);

end

